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Abstract 

Exact-exchange self-consistent calculations of the Kohn-Sham potential, surface energy, and work 
function of jellium slabs are reported in the framework of the Optimized Effective Potential (OEP) 
scheme of Density Functional Theory. In the vacuum side of the jellium surface and at a distance 
z that is larger than the slab thickness, the exchange-only Kohn-Sham potential is found to be 
image-like (~ —e^/z) but with a coefficient that differs from that of the classical image potential 
yim{z) = — e^/4z. The three OEP contributions to the surface energy (kinetic, electrostatic, and 
exchange) are found to oscillate as a function of the slab thickness, as occurs in the case of the 
corresponding calculations based on the use of single-particle orbitals and energies obtained in the 
Local Density Approximation (LDA). The OEP work function presents large quantum size effects 
that are absent in the LDA and which reflect the intrinsic derivative discontinuity of the exact 
Kohn-Sham potential. 



1 



I. INTRODUCTION 



The analysis of the electronic structure of metal surfaces poses a big theoretical challenge: 
a suitable calculational tool is needed for large, interacting, and strongly inhomogeneous 
many-electron systems. More than thirty years since its first application by Lang and Kohn 
to the surface problem,-!^ little doubt exists that one method of choice for the fulfilling of 
this goal is Density Functional Theory (DFT).^'^ DFT aims to a microscopic understanding 
of atoms, molecules, clusters, surfaces, and bulk solids starting from the fundamental laws 
of quantum mechanics. In the Kohn-Sham (KS) implementation of DFT,- the complicated 
many-body problem is mapped to an effective single-particle problem, with particles sub- 
jected to an effective single-particle potential (the KS potential). Although this mapping is 
exact, it gives no clue as to how to calculate in practice the so-called exchange-correlation 
(xc) contribution to the KS potential. Lang and Kohn solved this problem by using the 
Local Density Approximation (LDA) for the surface problem.-i^ In LDA, the xc potential 
at each point is taken to be that of a homogeneous interacting electron gas with the local 
density. Since then, many authors have calculated the electronic properties of metal sur- 
faces by using either the LDA~ or further elaborations that incorporate non-local ingredients 
to the unknown xc functional.-!^ Other schemes of the computational electronic-structure 
tool kit available for the investigation of solid surfaces are the Fermi hypernetted chain 
(FHNC) method,-SiiS the GW approximation,^ Quantum Monte Carlo (QMC),-^^'^^'^^ and 
the inhomogeneous Singwi-Tosi-Land-Sjolander (ISTLS) approach.— 

In the framework of the Optimized Effective Potential (OEP) scheme of DFT,— which 
had been first used in the context of atomic physics,— correlation is ignored altogether and 
the exact-exchange KS potential is obtained. Several advantages are associated with the use 
of the exact-exchange energy functional of DFT: (i) it corrects the self-interaction problem 
inherent in approximate treatments of the exchange energyi^ (this problem is particularly 
acute for localized systems such as atoms and molecules, although it is not relevant for 
extended systems like bulk solids and solid surfaces); (ii) it yields great improvements in the 
study of the KS eigenvalue spectrum,— semiconductor band structures and excitations,— 
and nonlinear optical properties;— (iii) it yields the correct asymptotics;— (iv) it reproduces 
the derivative discontinuity which should be present in the KS exchange potential each time 
the number of particles crosses through an integer value;— >2^'2^>2Ii^ and (v) it yields the 
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correct two-dimensional (2D) exchange energy per particle in the case of a quasi-2D electron 
gas.— It is the aim of this paper to provide benchmark exact-exchange OEP calculations for 
jellium slabs, with the expectation that more accurate DFT schemes that include correlation 
be developed by starting from a well founded exchange analysis and tested once reduced to 
their exchange-only (x-only) counterparts. 

The rest of the paper is organized as follows: We give in Section II the general theoretical 
background which will be used in the following sections; Section III is devoted to a discussion 
of the asymptotic behaviour of the exact-exchange KS potential of jellium slabs; in Sections 
IV and V we give the results that we have obtained for the OEP surface energy and work 
function, respectively, and in Section VI we present the conclusions. 

II. THE OEP APPROACH 

Our calculations are restricted to a jellium-slab model of metal surfaces, where the discrete 
character of the positive ions inside the metal is replaced by a uniform distribution of positive 
charge (the jellium). The positive jellium density is defined as 



which describes a slab of width d, number density n,~ and jellium edges at z = —d and 
z = 0; 9{x) represents the Heaviside step function: 6{x) = 1 if a; > and 6{x) = if a; < 0. 
A schematic view of our jellium slab is given in Fig. 1. Besides, and for convenience for the 
numerical calculations, infinite barriers are located far from the jellium edges, well inside 
the left and right evanescent vacuum regions. We have checked that these infinite barriers 
are located far enough for all the numerical calculations presented here to be independent of 
their precise location. The jellium-slab model is invariant under translations in the x — y 
plane, so the KS eigenfunctions can be factorized as follows 



where p and k are the in-plane coordinate and wave- vector, respectively, and A represents 
a normalization area. C,i{z) are the normalized spin-degenerate eigenfunctions for electrons 
in slab discrete levels (SDL) i {i = 1,2,...) with energy Si. They are the solutions of the 
effective one-dimensional KS equation 




(1) 




(2) 
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+ Vks (z) - Si 



Uz) = 0, 



2me dz'^ 

with rrie the bare electron mass. 

The KS potential Vks entering Eq. ([3]) is the sum of two distinct contributions: 

VKs{z) = Vn{z) + V,,{z), 
where Vh(-z) is the classical (electrostatic) Hartree potential, given by^ 

/oo 
dz'\z-z'\ [n{z') - n+{z')] . 
-00 

Here, n[z) is the electron number density— 



(3) 



(4) 



(5) 



(6) 



where kp = A/2me(/i — ei)/h, and /i = jj,{n,d) is the chemical potential, which in turn is 
determined from the neutrality condition for the whole system by the condition YlT'^' {^f)"^ ~ 
27r(in. Vxc{z) is the nonclassical xc potential, which is obtained as the functional derivative 
of the so-called xc energy functional E.ji.(.[n{z)\M' 

1 5E,,[n{z)] 



V.Az) 



(7) 



A 6n{z) 

Applications of DFT typically proceed from explicit density-dependent forms of E^c, 
as obtained using a variety of local or semi-local approximations. However, in the last 
few years increasing attention has been devoted to orbital-dependent forms of E^c'- E^c = 
Exc{{ii} 1 {siY\i which are only implicit functionals of the electron density n{z). In this 
case, one resorts to the OEP method^ or, equivalently, uses repeatedly the chain rule for 
functional derivatives to obtain the following expression for the xc potential of Eq. ([7]):— 



V^z) 



^ occ. 



dz' / dz 



[6Uz")6VKsiz') 



c.c. 



SVks{z') 
6n{z) 



(8) 



Multiplying Eq. (IHl) by the KS density-response function XKsi^^^') = ^''^{z)/6Vks{z'), using 
the identity 

00 

XKsi^^ z') Xks(^'' ^") dz = b{z- z), (9) 



comparing Eqs. ([7]) and ([8]), and integrating over the coordinate one finds 



6n{z' 



+ c.c. 



dz' . 



(10) 



The nice feature of Eq. lfTOl) is that 6C,i{z')/6VKs{z) and XKsi^^^') simply obtained 
from the solutions of Eq. ([3]) , as follows 



SVks{z) 



Uz) E 



Uz) Gf\z',z) 



(11) 



and 
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6n{z) 6^i{z") 
SUz")SVks{z': 



+ c.c. 



— oo 
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^^W*^^(^')^^'(^''^) + ^•^• 



(12) 
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where Gf^{z', z) is the Green function of noninteracting KS electrons. In the calculation of 
XKsi^i^')i chain rule for functional derivatives has been used; now we are considering 
the density itself as a functional of the occupied SDL. In obtaining Eq. f|T3|) from Eq. f|T2l) . 
we have used Eq. (ITT!) and also that 

5Uz') ^ ^ 27r ^'^^> ' 

which follows from Eq. ([6]). 

Introducing Eqs. ( fTTI) and ( fT3|) into the central Eq. ( JTOl) . we obtain the final and compact 
version of the OEP integral equation for Vxc{z): 



where 



and 



J2s,{z)=0, 

i 

Si{z) = {k'FY^.,{zyuz) + c.c. 

oo 

Uz) 



(14) 



UzrAv:,{z')Uz')dz'. 



(15) 



Here, AV^^{z) = Vxc{z)—ul.^{z), where ul.^{z) are SDL-dependent xc potentials of the form:— 

= [^7r/A{k'p)%{zr] 5ExJ5i, 
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The magnitudes "^i^z) are called the "shifts", as they can be physically interpreted as 
the first-order corrections of the KS eigenf unctions C,i{z) under the perturbation AV^^{z). 
These shifts also provide a useful and practical tool for the numerical solution of the OEP 
equation.— 1^ From Eq. f[T^ . we find the orthogonality constraint between the KS eigenfunc- 
tions and the shifts: J ^i{z)* "^{{z) dz = 0. It is also immediate that the shifts are invariant 
under the replacement Vxc{z) Vxc{z) + a, with a being an arbitrary constant. This means 
that the above set of equations determines Vxc{z) up to an additive constant, which should 
be fixed by imposing a suitable boundary condition. Moreover, the shifts "^iiz) are easily 
found to satisfy the following inhomogeneous differential equation:— 



AVUz) - AV: 



Uz) 



(16) 



Here, mean values are defined as O = / ^i{z)*0\z) ^i{z) dz. 

Equations ([3])- ([5]) and f|T^ . which determine the local Vxc{z) corresponding to a given 
SDL-dependent E^c, form a closed system of equations (the OEP equations), which should 
be solved in a self-consistent way. In order to accomplish some contact with other useful 
versions of the OEP equations for the present problem, a few additional steps are required. 
First of all, we write 



2mp 



Uzy 



d^^HAz) 



d%{z) 



dz'^ dz'^ 

which is easily obtained from Eq. Q. Secondly, we multiply the left hand-side of Eq. ([16 
by ii{z)* to obtain 



Uz) 



^^^^{Z) 

dz"^ 



"^iiz 



d%{zr 

dz"^ 



/Wliz) - AV' 



Uz)f 



Then, we start from the self-evident identity 



Vxc{z) 



47m{z) 



uUz) + AVl + AVUz) - AVl + c.c. 



(17) 



(18) 



we eliminate the factor \AV^^{z) — AV^^ 



\^i{z)f by using Eq. (fTTIl . and we obtain 



VUz) 



ATcn(z) 



Uz) 



uUz)+AV: 



+ 



2mp 



Uz) 



dz"^ 



^dz 



dz"^ 



+ c.c. 



(19) 



Finally, we proceed with the elimination from Eq. (fT9l) of the term proportional to 
d'^'^i{z) / dz^ , the subsequent elimination of d'^^i{z)* /dz^ proceeds via the KS equations. 



and as a result of all these manipulations we obtain the following expression for the DFT 
xc potential:— 

V.c{z) = V,,^i{z) + V,,^2{z), (20) 



where 



and 



Ke,2(^) = - 5^ (/i - [(fci.)' Vl/,,(^)e,(;^)* + VI/:(Z)C'(Z)* + C.C. 

with primes denoting derivatives with respect to the z coordinate. It is important to note 
that Eqs. f|T^ and fl20|) are just two different, but fully equivalent, ways to obtain the 
OEP xc potential for the present problem. If the shifts "^iiz) are (arbitrarily) forced to be 
identically equal to zero, the only term that survives is Vxc,i{z). This is exactly the KLI 
approximation,— which brings the identification Vxc,i{z) = V^^^z)."^"^ As before, Eqs. ([3])- 
([5]) and ( l20l) form a closed set of equations, which should be solved self-consistently. 

Both exchange and correlation have been included so far. Unless stated otherwise, we 
will now focus on the x-only case, where Er^c, Krc(^), and u^dz) are replaced by E^, V^^z), 
and Ux{z), respectively. We have achieved the self-consistent numerical solution of the x- 
only version of the OEP equations by two different methods: i) direct calculation of the 
shifts of Eq. f|T5l) . by solving Eq. f|T6l) .— and ii) direct solution of the OEP integral equation 
for Vx{z), as given by the a;-only version of Eq. (|T^ .'^° Both methods yield results that 
agree within numerical accuracy, although the first approach is found to be computationally 
more efficient than the second. Both methods face numerical instabilities beyond a critical 
coordinate z in the vacuum region. 

Finally, we note that the exact-exchange energy of a jellium slab is given by the following 
expression: 

oo 

^^(^) = i^E (^^)' J dzUz)\'ul{z), (21) 

-oo 

where ul.{z) represent the SDL-dependent exchange potentials 

oo 

^e,(^)* f ^ ,C.{zr9{^zk^,Azyp)^,{z') 
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with 



\z — z 



oo 

g{s, s') = j 



Ji{st)Ji{s't) dt 



(23) 





being the "universal" (that is, independent of Vks) function introduced by Kohn and 
Mattsson,^^ and Ji{x) being the first-order cyhndrical Bessel function.— 



III. ASYMPTOTICS OF THE EXACT-EXCHANGE KS POTENTIAL 



The long-range behavior of V^ciz) in the vacuum region is an important and open issue in 
DFT studies of metal surfaces.— The aim of this section is to present a detailed derivation of 
the analytical asymptotic limit of Vx{z) reported in Ref. |4ll for a slab geometry. First of all, 
we note that by making the choice that Vks{z —>■ oo) — > 0, Eq. Q leads us to the conclusion 



that C,i{z oo) 



: V-2me Si 



/'^ for all occupied % (disregarding a factor involving powers 



of z\ We also remark the following points: i) Due to the exponential decay of Vh(;z oo), 
the assumption Vks (-2 oo) implies that V^^z oo) 0; ii) for this choice of the 
zero of energy, one finds < for all occupied states; the slowest decaying of all the 
occupied SDL corresponds to i = m, where m is the highest occupied SDL. 

Now we look at the asymptotic behavior of the shifts "^{{z). Turning to the a;-only version 
of Eq. ([ISD, 



+ Vn(z) + VJz) - Si 



^,iz) = -VxizM^) + uliz)Uz) + AVM^), (24) 



2me dz"^ 

we focus on the asymptotic behavior of the three terms on the r.h.s. of this equation:— 

Vx{z oo) ii{z ^ oo) ^ 14(2; ^ 00) e"^^' 



ul{z 00) S,i{z 00) 
AF:e.(^ - oo)^e-^^% 



-z/3„ 



(25) 
(26) 
(27) 



with (3i = ^J—2mf,ei/%. Eq. fl26|) follows from an inspection of Eq. fl22l) in the limit z 00: 
in this limit, the sum over j is exponentially dominated by the term j = m, and the result 
of Eq. ( l26l) follows at once. Hence, for i ^ m Eq. ( 12^ yields 



2me dz"^ 



- Si 



"^Az 



00 



(28) 
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i.e., — > oo) e-^/3m^ Pqj^ ^ _ g^jj three terms in the r.h.s. of Eq. decay equally 
(to exponential accuracy), and further analysis is necessary. Eq. f|T4|) can be rewritten as 
follows 

m—l 

{k^Y ^mi^rUiz) + c.c. = ~Y1 (^f)' W*6 W - C.C., (29) 



1=1 



and by studying its asymptotic limit, it is clear that its r.h.s. can be approximated by the 
term i = m — l (with exponential accuracy). Given that both C,m{z oo) and \E'm-i(-2 oo) 
decay as e~^^"^'\ it follows that \l'm(-2 — ^ oo) decays as ^rn-i{z oo), that is, \£'m(-2 oo) — >• 
g-2;/3„_i^ Armed with these results, the asymptotic limit of Vx{z) is immediate from Eq. (1201) : 
Vx,2{z oo) tends exponentially to zero, while 

K,i(2 ^ oo) ^ u'^iz ^ oo) + ATT. (30) 

The leading contribution to u^{z oo) is easily obtained from Eq. fl2^ . by considering 
once again that in this regime the sum over j is exponentially dominated by the term j = m. 
For this case, the integral over the coordinate t can be evaluated analytically, yielding 



2 



M^(z ^ oo) ^ -e^ 



|^m(^ / I t / 

az 



\z — z 





z-z'\) 


Km 


Z 


- z' 





k'^\z~z'\ 



(31) 



-oo 



where Ii and Li are the modified Bessel and Struve functions, respectively.— Noting now 
that in this regime k]}\z — z'\ ~ A;™ 2; 3> 1^^^ it is permissible to expand the integrand of 



Eq. (|3T|1 as follows 

u^{z 00) ^ j 



00 

^2 



I' dz' 



z' 2 _ , 1 

^ p ' 



(32) 



Using the normalization of the orbitals ^ra{,z)^ we obtain 



<(^^oo)^-^(l + ^ + ...), (33) 

with /3 = IJ^ — 2/ (tt/c^) .— Since the exchange potential Vx{z) has been chosen to vanish 
at large distances from the surface into the vacuum [V^(oo) = 0], Eq. (l30l) leads us to the 
important constraint 

AFr=Fr-<=o, (34) 

which fixes the undetermined constant in Vx{z) discussed above. All numerical results 
presented here have been obtained by using this constraint. From Eqs. (I30l) . (|33|) . and (l3l|) . 



we conclude that 

V^{z ^ oo) ^ ^ oo) ^ ^ oo) ^ (^1 + ^ + ...^ , (35) 

which is the main result of this Section. 

At this point, we emphasize that the asymptotics dictated by Eq. (l35il hold only at z 
coordinates that are larger than l/k^. As A;™ is of the order oi 1/d (or smaller, depending 
on the actual value of rf), Eq. ( l35l) shows that the x-only KS potential happens to be four 
times larger than the classical image potential {Vimiz) = —e^/Az) only at a distance z 
that is considerably larger than the slab thickness. Furthermore, the arguments leading to 
Eqs. ( !26ll and (|3T1) are only valid for a discrete slab spectrum, such that there is a finite 
energy gap between Em and the remaining occupied energy levels Ei {i < m) . An extension 
of the present OEP framework to treat the case of a semi-infinite jellium surface^ is now in 
progress^. 

Finally, we note that under the condition l^s(oo) = Eq. fl35l) for the asymptotics of 
Vx{z) remains valid when correlation is included in the evaluation of the shifts "^{{z). The 
point here is that the shifts are separable in their exchange and correlation components, and 
they also satisfy separated differential equations (like Eq. (1241) for exchange). Once exchange 
and correlation contributions are splited, the analysis of the asymptotic behavior of Vx{z) 
follows the same lines as above, and the asymptotic limit of Eq. fl35|) remains the same. 



IV. SURFACE ENERGY 



In this section, surface-energy calculations are presented, as obtained at the x-only level. 
The surface energy a is the work required, per unit area of the new surface formed, to split 
the crystal in two along a plane.- For our slab geometry, 

^ '-m_im, (36) 

where E{d) is the total ground-state energy for each half of the slab after it is split (width 
d), and E{2d) is the total ground-state energy of the unsplit slab (width 2d), both the split 
and unsplit systems with the same jellium density. 

Following the standard DFT energy-functional partitioning, the surface energy (without 
correlation contribution) can be written as the sum of three terms,- 

a{d) = oxid) + aei{d) + ax{d), (37) 
10 



where (JK{d) is the non-interacting kinetic contribution to the surface energy, (Jei{d) is the 
electrostatic surface energy due to all non-compensated positive and negative charge dis- 
tributions in the slab, and (Jx{d) is the exchange contribution to the surface energy. From 
elementary physical arguments, it follows that axid) < 0, while (Jei{d) and cr^((i) are both 
positive.^ Also, the stability of the slab against spontaneous fragmentation is accomplished 
if cr(c/) > 0. From Eqs. (!36|) and (!37|) . one writes 



aiid) = [2Ei{d)-Eii2d)]/i2A), 



(38) 



with / = e/, X, and 



4 fc2 



Eei{d) = ^ I Vniz) [n{z) - n+{z)] dz, 



(39) 



(40) 



and [see Eqs. (l2T|-(j22l)] 



oo oo 



EM 



2n 



dz dz' 



(Az) 



(41) 



-oo — oo 



The dependence on the slab width d in Eqs. fl39|) . fl40l) . and fHTj) enters through the self- 
consistent KS eigenvalues {ei) and eigenfunctions {C,i{z)). 

Alternatively, one can define the effective single-slab surface energies^ 



a{d) 



E{d) - E'^'^'f {d) 



2A 



(42) 



and 



ai{d) 



EM)~Er'f{d)\/{2A), (43) 

where E^^^^ [d) is the ground-state energy of a uniform slab of electron gas of size d, 
and / = e/, xP^ Notice that Eq. fH2|) only reproduces the surface-energy definition of 
Eq. fl36p as d ^ oo. However, for a correct extrapolation of finite-slab calculations to the 
infinite-width limit,— here we calculate numerically the three components of the surface 
energy from the single-slab Eq. fl43p . We have checked that the differences between surface 
energies obtained from Eq. (!36|) and Eq. (H2l) are quite small even for the narrowest slabs 
studied, and that both agree in the extrapolation towards the semi-infinite limit. 
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Being the ground-state density the basic ingredient of DFT, we found interesting to 
compare the differences between the different density profiles that we have obtained. We 
exhibit in Figure 2 the self-consistent electron density profiles that we have obtained within 
the a;-only LDA and OEP schemes for = 2.07 and d = SXp.^"^ It is expected that the 
amplitude of the difference between both densities diminishes as z approaches the slab 
center, where both n^^^{z — d / 2) and rf'^^[z ^ — d / 2) should approach n as d ^ oo. 
Fig. 2 shows that there are noticeable differences between both densities: n^^^{z) extends 
further into the vacuum region than which is a result of the LDA orbitals being 

more extended or "diffuse" than their OEP counterparts, and the amplitude of the Friedel 
oscillations near the surface is larger for n^^^{z) than for n^^^{z). We have found the same 
behaviour for other values of r^. 

Figure 3 shows the results that we have obtained for the slab kinetic surface energy, as 
a function of the slab width d, for = 2.07. As in the case of the electron density, we 
have performed these calculations within the x-only LDA and OEP schemes. In the LDA, 
the kinetic surface energy axid) (LDA) is obtained by introducing the x-only self-consistent 
LDA eigenfunctions i\^^{z) and eigenvalues e\^^ into the formally exact Eq. ( l39l) . In the 
OEP, the kinetic surface energy axid) (OEP) is obtained by using the same Eq. fl39l) but 
with the LDA eigenfunctions and eigenvalues replaced by their x-only OEP counterparts 
^^^^{z) and ef^^. The strong oscillations in both axid) (LDA) and aK^d) (OEP) are 
the result of the sequential filling of empty slab discrete levels as d increases. Maxima in 
(Txid) correspond to the onset for the filling of a new slab discrete level. For this particular 
case, and following the extrapolation procedure of Ref. |50|, we have obtained the inf 
width extrapolated surface energies crx(LDA) = —4832 erg/cm^ (as reported in Ref. |5C 
and ai^(OEP) = - 4720 erg/cm^. 

Figure 4 displays the results that we have obtained for the electrostatic contribution to 
the surface energy, as a function of the slab width d and for = 2.07, again within the x-only 
LDA and OEP schemes. The electrostatic surface energies (Jei{d) (LDA) and aei{d) (OEP) 
are obtained from Eq. ( l40i) by using either the x-only LDA electron density n^^^{z) or the 
x-only OEP electron density respectively. In this case, the onset for the filling of 

a new slab discrete level is always associated with a minimum. Following the extrapolation 



procedure of Ref. 



50 



we have obtained the infinite- width surface energies indicated by arrows 



in Fig. 4: crei(LDA) = 1172 erg/cm^ (as reported in Ref. l50l) and crei(OEP) = 1103 erg/cm 
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In Fig. 5, we show the results that we have obtained for the exact-exchange contribution 
to the slab surface energy, as a function of the slab width d and for = 2.07, again within 
the x-only LDA and OEP schemes. As in the case of the kinetic and electrostatic surface 
energies, exact-exchange surface energies (^^{d) (LDA) and (y.j;{d) (OEP) [both derived from 
the formally exact Eq. (1411) ] are obtained by using either the x-only self-consistent LDA 
eigenfunctions il'^^i^z) and eigenvalues e\^^ or their x-only OEP counterparts if^^{z) and 
eP^^. For comparison, we have also calculated standard LDA-exchange surface energies^ 

oo 

a^DA = ^ J dz n'^^'^iz) {^^'^[^^^^(z)] - ^^^(n)} , (44) 

— oo 

where £™*-^(r;,) is the exchange energy per particle of a uniform electron gas of density n: 
e""*-^(n) = —3e^{37r'^nY^^ / (An), and n^^^{z) represents the x-only LDA electron density. 

All cr.j.(d) (LDA), a^id) (OEP), and cr;^^^((i) exhibit the characteristic oscillatory be- 
haviour also shown by the other components of the surface energy. As in the case of the 
electrostatic surface energy, the onset for the filling of a new slab discrete level is asso- 
ciated with a minimum. Fig. 5 shows that while the LDA [see Eq. fj44l) ] considerably 
overestimates the exchange surface energy, which is a known result, the exact-exchange 
surface energy is not very sensitive to the actual shape of the single-particle orbitals and 
energies, i.e., to whether LDA or OEP orbitals are used. Following the extrapolation proce- 
dure of Ref. [Sy, we have obtained the infinite-width surface energies indicated by arrows in 
Fig. 5: al;^^ = 2767erg/cm2, a^{LI)A) = 2390 erg/cm^ (both as reported in Ref. [sO), and 
a^(OEP) = 2316erg/cm2. 

We have also computed kinetic, electrostatic, and exchange surface energies for other 
values of the electron-density parameter r^, and we have obtained the infinite-width extrap- 
olated results shown in Table I. A comparison of the LDA and OEP calculations presented 
in Table I shows that (i) LDA orbitals being more delocalized than the more realistic OEP 
orbitals, surface energies that are based on the use of LDA orbitals are too large relative 
to those obtained with the use of OEP orbitals, and (ii) the sum of kinetic, electrostatic, 
and exchange surface energies are not very sensitive to whether LDA or OEP is used in the 
evaluation of the single-particle KS eigenfunctions and eigenvalues 
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TABLE I: Infinite- width extrapolated results for exchange-only kinetic [ax (LDA) and ctk (OEP)], 
electrostatic [a^i (LDA) and a^i (OEP)], and exchange [o";^;'^'^, ax (LDA), and ax (OEP)] surface 
energies for different values or r^. a (LDA) and a (OEP) represent the sum of the corresponding 
exchange-only kinetic, electrostatic, and exchange surface energies. Empty entries in a^i (OEP) 
for the two largest studied are due to the fact that the corresponding magnitudes are so small 
that it is not possible obtain a reliable extrapolated value. Units are erg/cm^. 





ax (LDA) 


aK (OEP) 


Tel (LDA) 


aei (OEP) 


LDA 


ax (LDA) 


ax (OEP) 


a (LDA) 


a (OEP) 


2.00 


- 5707 


- 5579 


1390 


1317 


3131 


2726 


2649 


-1591 


-1613 


2.07 


- 4832 


- 4720 


1172 


1103 


2767 


2390 


2316 


-1270 


-1301 


3.00 


-770 


- 733 


189 


177 


707 


568 


535 


-13 


-21 


4.00 


- 169 


- 155 


49 


48 


243 


180 


161 


60 


54 


5.00 


-46 


- 39 


19 




105 


71 


59 


44 




6.00 


- 13 


- 9 


9 




52 


32 


23 


28 





V. WORK FUNCTION 

The work function W is the minimum work that must be done to remove an electron 
from the metal at zero-temperature. In the context of DFT, the rigorous expression for the 
work function for a slab of thickness d is^ 

iy(rf) = l^Ks(oo)-/i, (45) 

where /i is the chemical potential. We note that as we are considering an electron system 
that is infinite in the x - y plane, electronic relaxation effects after removal of one electron 
are infinitesimal. For a slab geometry, the work function becomes size- dependent through 
the chemical potential fi{n,d). We are imposing the boundary condition yKs(oo) = 0; ac- 
cordingly, W{d) = — > 0. Besides, the only energy of the full KS spectrum which has 
a physical significance is precisely the energy of the highest occupied level, which can be 
identified with /i.^^ The work function for a slab with = 2.07 and d = AXp is shown 
schematically in Fig. 1. For this particular case, nine SDL are occupied and fi is between 
the nineth and tenth SDL. 

Now we focus on the slab-width dependence of the work function. Figure 6 shows the 
result of the x-only calculations that we have performed within LDA and OEP [ly'"^^ and 
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pyOEPj _ 2.07. The weakly oscillating x-only W^^^{d) is equivalent to the slab-width 

dependent work function reported by Schulte a long time ago.^^ As discussed by Schulte, 
the oscillations in W^^^{d) are the result of a combination of the shift of the bottom of 
the slab potential well and an effective film thickness shift, both effects suffering from an 
abrupt change each time the number of occupied SDL changes by one. The important point 
here, however, is the much stronger oscillations found in our W'~'^^{d) calculations, whose 
explanation is provided now with some detail. 

First of all, we note that, strictly speaking, the OEP work function W'~'^^{d) exhibits 
discontinuities of large size each time a new SDL becomes infinitesimally occupied. The 
first discontinuity in Fig. 6 appears at the 1 SDL 2 SDL transition (for d < Xf/2), the 
second discontinuity appears at the 2 SDL —>■ 3 SDL transition (for d Xp), and so on. In 
order of clarify the source of such a discontinuous behavior, we have plotted in Fig. 7 the 
OEP exchange potential V^(z) for slightly increasing values of the slab width d, around the 
6 SDL — s> 7 SDL transition. Each slab width d is characterized by a "filling factor" of the 
last occupied SDL, which is defined as follows 



Hence, — > 0+ (implying /i corresponds to an infinitesimally small filling of the 

last occupied SDL {i = m), while am 1~ corresponds to the threshold of occupancy 
of the next SDL {i = m + 1). The key point here is the dramatic change in Vx{z) when 
passing from the slab thickness corresponding to ag = 1^ to the infinitesimally thicker slab 
corresponding to aj = 0^ (~ 10"^). The remaining curves have been obtained for slab widths 
corresponding to the seventh SDL being progressively occupied: As ay increases from 0"*" to 
1~, Vx{z) approaches the form it had at ag = 1^, both in depth and asymptotic behavior, 
the only difference being a lateral shift of Vx{z) to the right that is simply due to the larger 
value of d. 

Secondly, we note that the potential barrier that forms at the interface, right after the 
jellium edge on the vacuum side of the surface, exhibits both V^^i^z) and Vx^2{.z) contributions 
[see Eq. (I20j . so the KLI approximation (which sets Vx,2{z) = 0) cannot be used for the 
analysis of the characteristic discontinuous behavior of the work function. In all cases in 
Fig. 7,Vx{z oo) 0. While this is clearly seen in the figure for the curves corresponding 
to ag = 1~ and = 1^ [in which case kp ~ l/d; see the asymptotics of Eq. fl35|) ]. it is 
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not evident at all for the set of potentials with small occupancies of the last occupied level, 
i.e., « 1. In this case, kp « 1/d and the asymptotic regime only takes place at z 
coordinates that go to infinity (as 07 0+) far beyond the z coordinates considered in 
Fig. 7. This is the situation for aj ^ 10^^, 10^^, and 10^^. As a final remark on this figure, 
it is important to realize that in the bulk and near the interface the exchange potentials 
Vx{z) corresponding to = 1~ and ay = 0"*" are simply related through a single vertical 
(constant) shift. This property, which can be verified numerically from Fig. 7, may also 
be derived analytically (see below). Finally, we note that although we have restricted our 
discussion to the case of a particular SDL transition, the same happens at every highest 
occupied lowest unoccupied SDL transition. 

With the aim of understanding how this discontinuous behavior of Vx{z) versus the slab 
width explains the results of Fig. 6 for the work function W'^^^{d), we show in Fig. 8 the 
slab OEP electronic structure just before occupation of the SDL ^ 7 (left panel), that is, 
at the slab width corresponding to 1~, and just after occupation of the SDL ^ 7 

(right panel), i.e., at the slab width corresponding to aj —>■ 0^. We note that while the 
Hartree potential approaches zero outside the surface exponentially and remains essentially 
unaffected by the infinitesimal population of the SDL #7 (compare left and right panels 
of Fig. 8), the OEP exchange potential (and therefore Vks{z) as well) suffers the abrupt 
jump explained in Fig. 7 which induces in turn the corresponding abrupt jump in the Fermi 
level. The net result in going from the left to the right panels of Fig. 8 is that the work 
function W^^^{d) suffers an abrupt (discontinuous) decrease, as the boundary condition 
l^s(oo) = is rigorously valid in both cases. This discontinuous behavior of W'^^^{d), 
shown schematically in Fig. 8, represents precisely the origin of the jumps that are visible 
in Fig. 6 at every threshold for SDL occupation. It is evident from Fig. 6 that the size of 
the discontinuity decreases as d increases. 

Finally, we investigate the size of the discontinuities that are visible in Fig. 6. For this, 
we rewrite the central OEP equation [as given by Eq. ([HI)] in the following way: 

m— 1 00 

Efe)' / mz';m)-uUz';m)]Gf''{z,z')v,{z',z)dz'+ 

00 (47) 
(fc-)2 / [V,{z'; m) - uliz'; m)] Gl^'iz, z') ^U^', z) dz' + c.c. = 0, 

—00 

where (pi{z, z') = C,i{z)*^i{z'). In writing Eq. (H71) the contribution of all the m — 1 occupied 
SDL's has been split from the contribution of the last occupied (m) SDL. The label m in 
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Vx{z; m) and m) has been introduced in order to emphasize that they are solutions of 
a system with m occupied SDL's. 

Let us now define a distance Z, such that ioi z > Z the electron density is dominated 
by the contribution of the last occupied (m) SDL, which is the one with the slowest decay. 
Eq. ([2]) clearly shows that Z ^ oo when k]^ 0, which is the case whenever — ^ O"*", 
i.e., whenever the filling of the last occupied SDL is infinitesimally small. We consider the 
following trial solution of Eq. fHTj) : 

V,{z;m) = V,iz;m-l) + Cx{m), (48) 

for z < Z and k]^ —* 0, with Cx{m) being a constant which depends on the last occupied 
SDL. Introducing this trial solution into Eq. (H7|) . we obtain 

m—l oo 

E [k^F? I m - 1) + Cx{m) - uliz'; m)] Gf^{z, z') ^,{z\ z) dz'+ 

oo (49) 
(^F)' / [VA^'; + - uiiz'; m)] G^z, z') ^^{z' , z) dz' + c.c. = 0. 

— oo 

In the limit —>■ 0, the second-line contribution of Eq. fj49l) is arbitrarily small; also, 
the KS wave-functions C,i{z) and eigenvalue differences (denominators) entering Gf^{z,z') 
should be extremely similar for the slab width corresponding to m — 1 occupied levels 
and am-i 1^, and the slab width corresponding to m occupied levels and am 0+. 
Therefore, an inspection of Eq. (!22|) leads us, using similar arguments, to the conclusion 
that u^^z^m) u^^z^m — 1), for all z < m, z < Z, and A;™ 0. Under these conditions, 
the first line of Eq. (H9l) reverts to the OEP equation for a slab width corresponding to 
m — l occupied states, and the proposal of Eq. ( l48l) is proved. Considering now that 
Cx{fn) = Vx{z;m) — Vx{z;m — 1), and taking the expectation value at the last occupied 
state (m — 1) of m — 1 system, we find 

Cx{m) = Vr\m) - Vr\rn - 1). (50) 

Now, for the m — l system we can use the boundary condition ^^^-'^(m — 1) = u^^^{m — l), 
and once again, approximate u'^~^{m — 1) ^ M™~^(m), yielding 

Cx{m) = Vr\m)-ur\m), (51) 

which has the nice feature that both the exchange potential Vx and the orbital-dependent 
exchange potential Ux are referred to the m system. For the m system V^{m) = u^{m), 
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TABLE II: Infinite-width extrapolated a;-only LDA and OEP work functions for various values of 
rg. Units are eV. 





2.00 


2.07 


3.00 


4.00 


5.00 


6.00 


p^LDA 


2.82 


2.80 


2.50 


2.15 


1.86 


1.62 


^OEP 


2.64 


2.63 


2.49 


2.11 


1.84 


1.61 



which does not prevent the constant Cx{m) from being nonzero (as shown in Fig. 7) since the 
KS orbitals C,m-i{z) and C,m{z) are different. As the slab width increases, m also increases 
and the difference between ^rn-i{z) and ^m{z) decreases, thereby leading to the expectation 
that Cx{m) — i> as (i ^ oo. This is explicitly shown in Fig. 9. While this analysis explains 
why Cxim) ^ for any finite m, it does not gives a hint about its sign; Fig. 9 shows, 
however, that Cx{m) is positive for all m. This positive jump in Vx{z) is exchange driven: 
at each threshold width for the occupation of a new level, a barrier appears against the 
occupancy of an empty SDL. This is due to the fact that intra-SDL exchange is stronger 
than inter-SDL exchange. As a consequence, the slab gains exchange energy by restricting 
new SDL occupancies. On the other hand, correlation induces in general a negative jump in 
Vc{z), so the net jump in Vxc{z) depends on the relative weigth of exchange and correlation 
for each particular system.— 

Finally, we have observed numerically that the average of the OEP work functions for 
slab widths corresponding to Om-i —^1^ and 0"*" remains the same (within error 

bars) for all the m values that we have considered. Hence, we have taken the infinite-width 
extrapolated work function to be simply that average. Table II exhibits the infinite-width 
x-only LDA and OEP work functions that we have obtained in this way for various values 
of the electron-density parameter r^. OEP work functions are slightly and sistematically 
smaller than their LDA counterparts. 



VI. CONCLUSIONS 

We have reported benchmark exact-exchange self-consistent calculations of the KS poten- 
tial, surface energy, and work function of jellium slabs in the framework of the OEP scheme. 
Special emphasis has been put into the asymptotical behaviour of the exact-exchange KS 
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potential far into the vacuum and the large quantum size effects that are present in the 
slab-width dependence of the surface energy and work function. 

We have performed a detailed analysis of the asymptotics of the exact-exchange KS 
potential far into the vacuum^, showing that at a distance z that is larger than the slab 
thickness the exact-exchange potential takes an image-like form: V^i^z — > oo) —e^ / but 
with a coefficient that differs from that of the classical image potential Vim{z) = —e'^/Az. 
Although this result has been obtained in the x-only approximation, it is also true in the 
presence of correlation due to the separability of the basic OEP equations in their basic 
exchange and correlation components. 

The OEP kinetic, electrostatic, and exchange contributions to the surface energy of jel- 
lium slabs have been obtained as a function of the slab width d and for a set of electron 
densities characterized by the parameter r^. We have shown that these components of the 
surface energy are all oscillating functions of d, with the oscillating period being ^ Ai?/2. By 
a suitable extrapolation procedure, we have found the values of the different components of 
the surface energy of a semi-infinite jellium. We have compared our OEP surface energies 
with those obtained from the same formally exact expressions [see Eqs. ( l39l) -( l4T|) ] but using 
single-particle LDA wave functions and energies; we have found small differences between 
these OEP and LDA surface energies, which appear as a consequence of the LDA orbitals 
being slightly more delocalized (diffuse) than their more realistic OEP counterparts. 

Finally, we have performed x-only OEP calculations of the work function of jellium slabs, 
again as a function of the slab width d. We have found that the OEP work function exhibits 
large quantum size effects that are absent in the LDA and which reflect the intrinsic deriva- 
tive discontinuity of the exact KS potential. The amplitude of this discontinuity diminishes 
as the slab width increases, and becomes arbitrarily small a.s d —>■ oo, i.e, in the case of 
a semi-inflnite system. This has been proved both analytically and numerically. We also 
note that although the precise value of the x-only OEP work functions reported here would 
change with the inclusion of correlation, the exact slab work function is expected to exhibit 
the large quantum size effects and discontinuities observed in the present work, barring pos- 
sible accidental cancellations of exchange-driven and correlation-driven contributions to the 
total discontinuity. The presence of these large discontinuities in the x-only OEP slab work 
function (and presumably also in the actual work function that includes correlation) high- 
lights the potential danger in which can be incurred by performing elaborated calculations 
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for a restricted set of slab sizes without performing a suitable and reliable extrapolation 
towards the semi-infinite case. 

In summary, we expect that the benchmark exact-exchange OEP calculations reported 
here for jellium slabs will serve as motivation and as a starting point for the development of 
more realistic approximations for the exchange-correlation energy functional of jellium and 
real surfaces. 
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FIGURE CAPTIONS. 



Figure 1. Main features of the jellium-slab model of metal surfaces. Top panel: 
normalized jellium density (n+(z)), and the self-consistent OEP electron density n{z) for 
two different values of the electron-density parameter Vg. Lower panel: OEP Hartree, 
exchange, and Kohn-Sham potentials for = 2.07. Dotted lines denotes KS eigenvalues, 
ep is the Fermi energy, and W is the work function, d = 4 Xp- 

Figure 2. LDA and OEP self-consistent electron densities and its difference for 
d = 8 Xp and = 2.07. Note that n^^^{z) is slightly more diffuse than n'^^^{z), as 
n^^^{z) — n°^^(z) > for z outside the jellium edge (in the vacuum). 

Figure 3. Kinetic surface energy, as a function of slab width d, for = 2.07, from 
Eq. ( I43l) . with / = K. Full line, OEP results; dotted line, LDA results. The two arrows 
on the right denote the extrapolated asymptotic values ax {OEP) —4720 erg/cm^, 
ax (LDA) ^ -4832 erg/cm^. 

Figure 4. Electrostatic surface energy, as a function of slab width d, for = 2.07, 
from Eq. (j^3|) . with / = el. Full line, OEP results; dotted line, LDA results. The two 
arrows on the right denote the extrapolated asymptotic values Uez (OEP) —>■ 1103 erg/cm^, 
aei (LDA) ^1172 erg/cm^. 

Figure 5. Exchange surface energy, as a function of slab width d, for = 2.07, from 
Eq. fH3|) . with / = x. Full line, OEP results; dotted line, LDA results; dash-dotted line, 
standard LDA-exchange results. The three arrows on the right denote the extrapolated 
asymptotic values a^(OEP) ^ 2316 erg/cm^, a^.(LDA) 2390 erg/cm^, a^^^ 2767 
erg/cm^. 

Figure 6. Slab work function versus slab width d, for = 2.07. Full line, OEP result; 
dotted line, LDA result. Occupation events corresponding to transitions from a slab with 
m occupied SDL towards m + 1 occupied SDL are denoted as m — > m + 1. 
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Figure 7. Self-consistent OEP exchange potential, around the 6 — 7 SDL transition, 
for Vs = 2.07. The origin of coordinate z for each slab has been taken at the slab center. 
The position of the right slab edge has been indicated by a vertical dashed line for each case. 

Figure 8. Left: electronic structure of the slab for ag = 1^. Right: electronic structure 
of the slab for aj = 0~^. The work function W jumps discontinously from its left large value 
towards the smaller right value. Slab edge is at z = 0. 

Figure 9. Exchange- driven discontinuity Cxim) for increasing number of occupied slab 
levels, as follows from Eq. EH 
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